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We report the results of a complete modal and nonmodal linear stability analysis of the 
electrohydrodynamic flow (EHD) for the problem of electroconvection in the strong injec¬ 
tion region. Convective cells are formed by Coulomb force in an insulating liquid residing 
between two plane electrodes subject to unipolar injection. Besides pure electroconvec¬ 
tion, we also consider the case where a cross-flow is present, generated by a streamwise 
pressure gradient, in the form of a laminar Poiseuille flow. The effect of charge diffu¬ 
sion, often neglected in previous linear stability analyses, is included in the present study 
and a transient growth analysis, rarely considered in EHD, is carried out. In the case 
without cross-flow, a non-zero charge diffusion leads to a lower linear stability threshold 
and thus to a more unstable flow. The transient growth, though enhanced by increas¬ 
ing charge diffusion, remains small and hence cannot fully account for the discrepancy 
of the linear stability threshold between theoretical and experimental results. When a 
cross-flow is present, increasing the strength of the electric field in the high-i?e Poiseuille 
flow yields a more unstable flow in both modal and nonmodal stability analyses. Even 
though the energy analysis and the input-output analysis both indicate that the energy 
growth directly related to the electric field is small, the electric effect enhances the lift-up 
mechanism. The symmetry of channel flow with respect to the centerline is broken due to 
the additional electric field acting in the wall-normal direction. As a result, the centers of 
the streamwise rolls are shifted towards the injector electrode, and the optimal spanwise 
wavenumber achieving maximum transient energy growth increases with the strength of 
the electric field. 


1. Introduction 

1.1. General description of EHD flow 

Electrohydrodynamics (EHD) is concerned with the interaction between an electric field 
and a flow field. Such configurations have broad applications in a range of industri al and 

biological devices. EHD effects c an be used to enhance the heat transfer efficiency ([Jones _ 

Allen fc KaraviannisIlQQ^ , to design microscale electrohydrodynamic pumps (iBart etHil. 


1978; 


199C; 


Darabi et a?.ll20021 . to fabricate diagnostic devices and drug delivery systems (jChakrabortv et al. 
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200it) a nd DNA microarravs dLee et Qi.ll2006h , and to design new strategies for active flow 
control ( Bushnell fc McGinle^ 1989 1. Physically, EHD flow is characterized by a strong 
nonlinear interaction between the velocity field, the electric field and space charges: the 
electric force results in flow motion, which in turn affects the charge transport. The 
intricate nature of this nonlinearity defies a fundamental understanding of EHD flow. 
Moreover, as we will see, there still remains a mismatch or discrepancy between experi¬ 
mental observations and a theoretical analysis. 

One classic problem in EHD, named electroconvection, deals with the convective mo¬ 
tions induced by unipolar charge injection into a dielectric liquid (of very low conduct¬ 
ivity) which fills the gap between two parallel rigid plane electrodes. The Coulomb force 
acting on the free charge carriers tends to destabilize the system. Electroconvection 
is often compared to Rayleigh-Benard convection (RBC) because of their similar geo¬ 
metry and convection patterns. Moreover, RBC is known to be analogous to the Taylor- 
Couette (TC) flow in the gap between two concentric rotating cylinders, where thermal 
energy transport in RBC corresponds to the transport of angular momentum in TC flow 
I Bradshaw Il969l : ICrossmann fc Lohsel 2000ll . In the linear regime of RBC, the flow is 
destabilized by the buoyancy force caused by the continued heating of the lower wall (an 
analogous role is played by centrifugal force in TC flow). As the thermal gradient ex¬ 
ceeds a critical value, chaotic motion sets in. In EHD flow, the destabilizing factor is the 
electric force, acting in the wall-normal direction. However, the analogy between the two 
flows ends, as soon as nonlinearities arise, especially as diffusive effects are concerned: 
in RBC, molecular diffusion constitutes the principal dissipative mechanism whereas in 
EHD flow, it is the ion drift velocity KF (with K being the ionic mobility) which diffuses 
perturbations in the fluid. It is well-known that RBC is of a supercritical nature, i.e., 
transition from the hydrostatic state to the finite-amplitude state occurs continuously as 
the controlling parameter, i.e., the Rayleigh number, is increased. For EHD flow, on the 
other hand, the bifurcation is subcritical, characterized (i) by an abrupt jump in motion 
amplitnde from zero to a finite value, as a critical parameter is crossed, and (ii) by the 
existence of a hysteresis loop. It is interesting to mention an analogy between EHD flow 
and polymeric flow: polymeric flow shows a hysteresis loop as well, as the first bifurcation 
is considered. In fact, the counterpart of EHD flow, i.e., magneto-hydro dynamics (MHD) 
flow, has been compared to polymeric flow in lOgilvie fc Proctor ( 2003h . 

Most studies in the EHD literature address electroconvection in the hydrostatic con¬ 
dition, i.e., without cross-flow. In this work we also investigate the EHD stability prop¬ 
erties in the presence of cross-flow. Our interest is two-fold. First, the potential of this 
flow configuration resides in the possibility of using the electric field to create large- 
scale rollers for flow man ipulation; turbulent drag reduction designed i n the spirit of 
SchoDpa fc Hussaii] ( 1998l l and investigated by Soldati fc Baneri^ (1998) in the nonlin¬ 


ear regime is an example of this type. Secondly, EHD with cross-flow has been applied to 
wire-plate electrostatic precipitators, but due to the complex nature of the chaotic inter¬ 
action between wall turbulence and the electric field, our current understanding of such 
flows is rather li mited. Nonlinear EHD sim ulations with a cross-flow component have 
been reported in Soldati fc Baneri^ ( 1998l l. More relevant to our linea r problem is the 
unipo lar-i nj ection-induced instabilities in plane parallel flows studied by lAtten fc Honda 
( 1982h and Castellanos fc AgraitI ( 1992 1. The former work focused on so-called electrovis- 
cous effects, defined by an increase of viscosity due to the applied electric field compared 
to the canonical channel flow. The latter work found that, at high Reynolds numbers, the 
destabilizing mechanism is linked to inertia, while, at sufficiently low Reynolds numbers, 
EHD instability are dominant. In this article, we will not only address the modal stability 
problem of EHD channel flow, as those two previous studies did, we will also take into 











































3 


Modal and nonmodal stability analysis of electrohydrodynamic flow 

account transient effects, discussed shortly below, of the high-i?e number channel flow 
in the presence of an electric field. Our results would be interesting to the researchers in 
the flow instability and transition to turbulence, especially for high-i?e flow. The results 
will also shed light on the study of flow control strategy using EHD effects. 


1.2. Stability of EHD flow 

The endeavor to underst and the stability and transi tion to turbulence in EHD flo w dates 
back to the 1970’s, when Schneider fc WatsonI ( 1970f) and Atten &: Moreau ( 1972li . among 
the first, performed a linear stability analysis on the flow of dielectric liquids confined 
between two parallel electrodes with unipolar injection of charges. The mechanism for 
linear instability could be explained via the formation of an electric torque engendered 
by the convective motion when the dri ving electric force is suf ficiently strong to overcome 
viscous diffusion. It was established in lAtt en fc More^ (jl972h that, in the weak injection 
limit, C <C 1, with C as the charge injection level, the flow is characterized by the criterion 
TcC^ ~ 220.7, where Tc is the linear stability criterion for the stability parameter T, 
defined in the mathematical modeling section [2j^ and, in the case of space- charge-limited 
(SCL) injec tion, C —» oo, they foun d ~ 160.75. However, according to lLacroix et al. 
1 19751) and Atten fc Lacroixl ( 1979f) . the experimentally determined stability criterion 
was notably diffe r ent fr om the theoretical calculations. In the experiments performed by 
Atten fc Lacroix 1 1979[) . the linear criterion was found to be Tc ~ 100 in the case of SCL, 


which is far lower than the theoretically predicted value. It was argu ed then that this 
disagreement might be due to neglecting charge diffusion ( Attenl[l976 b We will address 
this discrepancy in the SCL in this paper, and confirm that charge diffusion is indeed an 
important factor influencing the linear stability criterion in this c a se. 

The first nonlinear stability analysis was performed by F elicil (1971), who assumed 
a two-dimensional, a priori hydraulic model for the velocity field in the case of weak 
injection between two parallel plates. It was found that within the interval [T„;,Tc], 
where is the nonlinear stability criterion for T, two solutions exist, namely, a stable 
state and an unstable finite-amplitude state. This finding corroborated the fact that 
the bifurcation in the unipolar injection problem is of a su bcritical nature a n d tha t 
the flow has a hysteresis loop, as experimentally verified by Atten fc Lacroixl ( 1979f) . 
Physically, this subcritical b ifurcation is related to the formation of a region of zero charge 
( Perez fc Castellanoslll98^ . Later, this simple hydraulic mode l was extended to three - 
dimensional, hexagonal conyective cells for the case of SCL by lAtten fc Lacroixl ( 1979l i. 
and it was shown that the most unstable hydrodynamic mode consists of hexagonal cells 
with the interior liquids flowing towards the injec tor. The nonlinea r stabi lity criterion for 
three-dimensional, hexagonal cells, according to lAtten fc Lacrob] ( 1979I) . was Tni « 90 
in the experiments, but theoretical studies produced « 110. 

Most of the previous linear stability analyses of EHD flow focus on the most unstable 
mode of the linear system, which is insufficient for a comprehensive flow analysis. In 
fact, theoretically, the linear stability analysis is linked to the characteristics of the lin¬ 
earized Navier-Stokes (N-S) operator C which, in the case of shear flows (in this paper, 
the cross-flow case), may be highly nonnormal, i.e., C f- CC^ (with denoting 
the adjoint of C) or, expres sed differently, the eigenvectors of the linear opera tor are 
mutually nonorthogonal (see iTrefethen et al\ Il993 ; Schmid fc Henningsonll200lh . For a 
normal operator C = the dynamics of the perturbations is governed by the 

most unstable mode over the entire time horizon. In contrast, a nonnormal operator has 
the potential for large transient amplification of the disturbance energy in the early linear 
phase, ev en though the most uns t able mode is stable. The theory of nonmodal stability 
analysis ( Farrell fc loannoul 19961 : SchmidI 2007h . the main tool to be used in this work, 
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has been applied successfully to explain processes active during transition to turbulence 
in several shear flows. The fact that the bifurcation of EHD flow is subcritical, a trait 
often observed in shear flows governed by nonnormal linearized operators, tempts one to 
think that the discrepancy between the experimental value Tc ~ 100 and the theoret¬ 
ical value Tc ~ 161 in the SCL regime of EHD flow might be examined in the light of 
nonmodal stability theory. In fact, it seems surprising that this type of sta bility analysis 
has so far only rarely been applied to EHD flows, except for the work of AttenI 1 1974) 
in the case of hydrostatic flow. The method we employ here is different from Atten’s 
quasistationary approach: nonmodal stability theory, based on solving the initial-value 
problem, seeks the maximum disturbance energy growth over entire time horizon when 
considering all admissible initial conditions a nd identifies the optimal initial condition 
for achieving this maximum energy growth. In lAtteiil ( 1974h . a quasistationary approach 
was taken that proposed that disturbances grow rapidly, compared to the time variation 
of the thickness of the unipolar layer; however, transient energy growth due to the non¬ 
normality of the linearized operator in hydrostatic EHD has been found to be rather 
limited in this work. This is in contrast with EHD Poiseuille flow, where nonnormality 
is prevalent and should be considered from the outset. 

The present paper extends the work by iMartinelli et al. ( 2011 ) and is organized as 
follows. In § [21 we present the mathematical model, the governing equations and the 
framework of the linear stability analysis. In § [31 numerical details are given and a code 
validation is provided in the appendix. We then present in §|3]the results of the modal and 
nonmodal stability analysis and in §|5]the energy analysis. Finally, in §|6l we summarize 
our findings and conclude with a discussion. 


2. Problem formulation 

2.1. Mathematical modeling 

We consider the planar geometry sketched in figure [U where the Cartesian coordinate 
system used in this work is [x, y, z) or (l^,, ly, 1^) as the streamwise, wall-normal and 
spanwise directions, respectively. The two flat electrodes extend inhnitely in the x- and 
z-directions, and the applied voltage only varies in the y-direction. The distance between 
the two electrodes is 2L*. The dimensional variables and parameters are denoted with a 
superscript *. The electric field satisfies the reduced Maxwell equations. The c harges are 


gener ated through electrochemical reactions on the charge-injecting electrode (|Ali et al 


1985h . Since the electric conductivity is very low, conduction currents are negligible even 
in the presence of larg e electric fields. Therefore, mag netic effects in the Maxwell equa¬ 
tions can be neglected ( Melchetill981[ICastellanoslll998 1. leading to the quasi-electrostatic 
limit of the Maxwell equations 


dQ* 

dt* 


V* X E* 

= 0, 

(2.1a) 

V* D* 

= Q\ 

(2.16) 

-f V* • J* 

= 0, 

(2.1c) 


where E* is the electric field, D* = e*E* denotes the electric displacement, e* stands for 
the fluid permittivity which we assume constant here, Q* represents the charge density 
and J* is the current density. Considering equation (I2.1h l. it is a well-known practice 
to define a potential field (j>* according to E* = Combining the first two equa- 
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Figure 1. Sketch of the electroconvection problem with coordinate system (x, y, z). In the 
non-hydrostatic case, a flow rate is induced along the streamwise (x) direction. 


tions (j2.1h l and (I2.1b h we can write the governing eqnation for cj)* as 


7*2 1* 


' e* 


( 2 . 2 ) 


The current density J* arises from several sources. By modeling the EHD flow with only 
on e ionic species in a perfectly insulating fluid (conductivity a* =0), one can express J* 
as ( CastellanosI 1998h 

J* =K*E*Q*+U*Q*-DIV*Q* (2.3) 

where the first term accounts for the drift of ions (with respect to the fluid) under the 
effect of the electric field, moving at the relative velocity K*E*, with K* as the ionic 
mobility, the second term represents the convection of ions due to the fluid velocity 
U*, and the last term takes i nto account the charge diffu sion, with D* as the diffusion 
coefficient. Since the work of Eerez_fc_Castellanpa (1989l l. the vast body of literature, 
with the exception of Kourmatzis fc Shrimptonl (j2ni2l l to turbulent EHD flow, neglects 
the charge-diffusion term because of its small value when compared to the drift terms. 
However, we will show that, even though the numerical value of D* is very small, its 
impact on the flow dynamics is undeniable. 

The flow field is incompressible, viscous and Newtonian and governed by the Navier- 
Stokes equations, which, in vector notation, read 


,dU* 

dt* 


y* -u* = 0 , 

p*{u* ■ y*)u* = -v*p* + p*y*^u* 


(2.4a) 

(2.4&) 


where U* is the velocity field, P* the pressure, p* the density, p* = p*iz* the dynamic 
viscosity {ly* the kinematic viscosity) and E* the volumic density of electric force, which 
expresses the coupling between the fluid and the electric field. In general, E* can be 
written as 


1 


F* =Q*E* --\E*\^V*e* + V* 




dp*\ ’ 


(2.5) 


where the three terms on the right-hand side represent, respectively, the Coulomb force, 
the dielectric force and the electrostrictive force. The Coulomb force is commonly the 
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strongest force when a DC voltage is applied. As we assume an isothermal and homogen¬ 
eous fluid, the permittivity e is constant in space. As a result, the dielectric force is zero 
(however, it would be dominant in the case of an AC voltage). The electrostrictive force 
can be incorporated into the pressure term of the Navier-Stokes equation as it can be 
expressed as the gradient of a scalar field. Therefore, the only remaining term of interest 
in our formulation is the Coulomb force. 

The system is supplemented by suitable boundary conditions. In our problem, we 
assume periodic boundary conditions in the wall-parallel directions. The no-slip and 
no-penetration conditions for the velocities are assumed at the channel walls. For the 
potential held, we require Dirichlet conditions on both walls, on the injector (j}*{L*) = (/)q 
and the collector (jf{—L*) = 0 in order to hx the potential drop A(/)q between the 
electrodes. The injection mechanism is autonomous and homogeneous, meaning that the 
charge density is constant on the injector, not inhuenced by the nearby electric held 
and has a zero wall-normal hux of charge on the collector, i.e., Q*{L*) = —Qo 

= 0. Owing to the homogeneity in the wall-parallel directions, there is no 
requirement for boundary conditions in the x- and z-direction. 


2.2. Nondimensionalized governing equations 

In the no-crosshow case, as we are interested in the effect of the electric held on the 
how dynamics, we nondimensionalize the full system with the characteristics of the elec¬ 
tric held, i.e., L* (half distance between the electrodes), AcpQ (voltage difference applied 
to the electrodes) and Qo (injected charge density). Accordingly, the time t* is nondi- 
mensionalized by L*'^/{K*A4>q), the velocity U* by K* A(j)Q/ L*, the pressure P* by 
PgAr*^A(()g^/L*^, the electric held E* by Aejr^/L* and the electric density Q* by Qg. 
Therefore, the nondimensional equations read 


dt 


V • [7 = 0, 

^ + [U-W)U = -VP + 


dt 

V -[{E- 


U)Q] = 

V^(j} = -CQ, 
E = -V(/) 


CM^QE, 


where 


M = 




K* 


T = 




C = ^0^ 


A^oe* ’ 


Fe = 


K*A(j)l 

D* 


(2.6a) 

( 2 . 6 &) 

(2.6c) 

(2.6d) 

(2.6e) 


(2.7) 


Additionally, the nondimensional boundary conditions are I7(±l) = 0, = 1, <)>(—1) = 

0, Q(l) = -land^(-l) = 0. 

Various dimensionless groups appear in the equations as written above. M is the ra¬ 
tio between the hydrodynamic mobility (e/po)^ and the true ion mobility K. Gases 
us ually take on a value of M less than 0.1 and liquids have values of M greater than 
1 ( Castellanos fc AgraitlflOO^l . T (Taylor’s parameter) represents the ratio of the Cou¬ 
lomb force to the viscous force. It is the principal stability parameter, assuming a similar 
role as the Rayleigh number in Rayleigh-Benard convection. C measures the injection 
level. When C 3> 1, the system is in a strong-injection regime, and when C <C 1, it is in a 
weak-injection regime. Fe is the reciprocal of the charge diffusivity coefficient. The factor 
M'^/T appearing in equation (I2.6b i can be interpreted as the ratio between the charge 
relaxation time L*'^/{K*Acj)^) by drift and the momentum relaxation time jv*. This 
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mathematical model for EHD flow has been assumed and studied in many previous in¬ 
vestigations of the linear sta bility and turbulence analyses for a dielectric liquid subjec t 
to unipolar injection of ions ( Lacroix et al. I ll975t l^aore fc Perejl 201 ^IWu et'^l2013 L 
ex cept that the diffusion term in eq uation (I2.6b l is usually neglected (excluding the study 
of iKourmatzis fc Shrimotonl ( 2012n ). 


2.3. Linear stability problem 

The linear problem is obtained by decomposing the flow yariable as a sum of base state 
and perturbation, i.e., U = U + u, P = P+p, E = E + e, D = D + d, Q = Q + q 
and (j) = (j) + ip. For the vector fields, we have u = (u,v,w) and e = ( 61 , 62 , 63 ) along 
the three Cartesian coordinate directions. After substituting the decompositions into the 
governing equations (l 2 . 6 b -eL subtracting from them the governing equations for the base 
states and retaining the terms of first order, the linear system reads 


V • M = 0, 


(2.8a) 


du — — — — 

— + (m • V)U + (U ■ V)u = -Vp + —+ CM^{qE + Qe), (2.86) 


■[{E + U)q+{eP u)Q] = 


(2.8c) 

VV = -Cg, (2.8d) 

e = —V(/ 3 , ( 2 . 8 e) 

with the boundary conditions for the fluctuations m(±1 ) = 0 , </ 3 (±l) = 0 and q{l) = 
9q I 
dy ' 




2.3.1. Base states 

The base states are the solutions to equations (I2.6h -el in the case of no time depend¬ 
ence. Owing to the periodicity in the wall-parallel directions, we can reduce the shape 
of the base states as functions of y only, that is, U = U{y)ly and E = E{y)ly. For the 
base flow U{y), we are interested in the hydrostatic and pressure-driven Poiseuille flows 
which, after nondimensionalization, are given by 


U{y) = 0, 


U{y) = Re—{1 - y^) = (1 - y^), 


(2.9) 


respectively, in which the (electric) Reynolds number is defined as Re = = K* 

(in order to enforce the same constant flow rate). It is a passive parameter in the hy¬ 
drostatic case, but becomes a free parameter in the presence of high-i?e cross-flow, in 
which, consequently, M would be the passive parameter. Therefore, in the Poiseuille flow 
case, we modify the governing equation (I2.8b l by substituting the relation Re = T/M'^ 
to obtain 




1 


CT, 


{U ■ V)tt = — Vp -I- —-I- -—{qE + Qe). 


Re 


Re 


( 2 . 10 ) 


By doing so, it is more obvious to identify the effects of T and C on the electric force 
term. The parameter Re = K* Ac^q/v* here coincides with the canonical hydrodynamic 
equivalent Re^ = U*L*/v* because of the electric scaling we chose. However, in a general 
sense, the two may not necessarily be identical. The nondimensional quantity 

L*/U* 


Re K*A(l)l 

'r^ - u*L* ~ L*‘^/{K*A(j}*) 


( 2 . 11 ) 




















Figure 2. The base states: (a) U, (b) (/>, (c) Q. 


relates the eddy turn-over time and the charge relaxation time by the drift. According 
to the equality T = Re ■ , when Re is near linear criticality at 5772 and T is around 

10^, M w 0.1. It implies that the working liquid is gas. Moreover, in contrast to the 
nonlinear constitutive modeling for polymers in viscoelastic flow, the base flow is not 
modified under the influence of the base electric field, even though the coupling between 
U and Q is nonlinear in equation (j2.6b 'l. This is because the directions of the base flow 
and the base electric held are perpendicular. Nevertheless, the base pressure gradient in 
the wall-normal direction is no longer zero. 

The base electric held E{y) can be solved from equations (I2.6b -el. recast into an 
equation only for tp which reads 

0>'" + (^")^ + ^r' = O, (2.12) 

r e 

where prime ' denotes the spatial derivative with respect to the y-direction. The boundary 
conditions are ^(1) = 1, 0(—1) = 0, ^"(1) = —C and = 0. Analytical solutions 

to this fourth-order ordinary differential equation can be obtained by observing that 
the equation can be transformed into a Riccati equation; alternatively, as we do here, 
a simple numerical integration combined with a nonlinear gradient method provides us 
with the required (^( 2 /)-prohle. The Poiseuille base how and the base states of the electric 
and charge helds are shown in hgure [2] 

2.3.2. Matrix representation 

In linear stability analysis, it is a common practice to rewrite the huid system (I2.8h -bl 
in terms of the wall-normal velocity v and the wall-normal vorticity rj = dzU — d^w by 
eliminating the pressure term. For the electric held, the three equations (I2.8b -el can be 
reduced to one for (p. Therefore, the governing equations (j2.8h -el become, in terms of a 
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V — T] — ip formulation, 


dt 


drj 

avV 

dt 


- a , -d . 

+M2 [ - <^"'(V2 - ^)ip + 

-a -.dv M 2 2 

I I im I 1 yT4 


with boundary conditions 


u(±l) = 0,u'(±l) = 0, 
r 7 (±l) = 0, 

^(± 1 )= 0 ,^"( 1 ) = 0 ,^'"(- 1 ) = 0 . 


(2.13a) 

(2.13&) 

(2.13c) 


(2.14a) 

(2.146) 

(2.14c) 


For compactness, we write 7 = {v,ri,(f)'^, and the linearized system, recast in matrix 
notation, becomes 


/V 2 0 0 \ o /u\ /Los 0 Ly^\ /v\ 

0 I 0 - U = L, Ls, 0 U (2.15) 

Vo 0 vy \ifj 0 \ifj 

where I denotes the identity matrix and the submatrices Los, Ly^p, Lc, Lgq, L^y and 
Lpp can be easily deduced from equations (I2.13b ,-cl. To represent the system even more 
compactly, we can rewrite the linearized problem (12.1511 as 

= ^ ^ = <'■“> 

where C = A.~^B represents the linearized Navier-Stokes operator for EHD flow. 

Since the flow is homogeneous in the wall-parallel directions, the perturbations are 
assumed to take on a wave-like shape. Moreover, as we consider a linear problem with a 
steady base flow, it is legitimate to examine the frequency response of the linear system 
for each frequency individually. These two simplifications lead to 

f{x, y, z, t) = f{y, f) eyi-piiax -f ifz) = f{y) exp(-iwt) exp(zax -f if3z), (2.17) 


where / could represent any flow variable in {u,p,e,q,(p)’^, f{y,t) and f{y) are the 
shape functions, a and /3 are the real-valued streamwise and spanwise wavenumbers, 
and the complex-valued w is the circular frequency of the perturbation, with its real 
part LUy representing the phase speed and its imaginary part cui representing the growth 
rate of the linear perturbation. Upon substitution of the above expression into the linear 
problem (I2.16|l . we arrive at an eigenvalue problem for the v — y — ip formulation which 
reads 

— iuj = £ 7 , (2-18) 

where —iuj is the eigenvalue and 7 is the corresponding eigenvector. Both formulations, 
( 2 d^ and (|2.18l) , would be relevant as discussed in a recent review by I Schmid fc Branch 
( 2014^ . The least unstable eigenvalues obtained from the eigenproblem formulation (I2.18|) 
would determine the asymptotic behavior of the linear system, while the initial-value 
problem formulation (12.161) could be used to examine the dynamics of the fluid system 
evolving over a finite time scale. 
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2.3.3. Energy norm 

In our calculation of the nonmodal transient growth, we define the total energy density 
of the perturbation contained in a control volume O as 


J £*dV* = J ( 4 * + £;) dy* = J I {poU* ■u* + e*- d*) dV* 

= + w*^) + (2.19) 

The perturbed electric energy f * follows the definition in CastellanosI ( 1998ll . In terms of 
the v — rj — (^-formulation, the nondimensionalized energy norm in spectral space becomes 

^I -f -^D\Di 0 0 




= / 

JQ. 


0 

0 


idy 


M^{k^I + DlDi 


( 2 . 20 ) 


where the superscript ^ denotes the complex conjugate, -1-/3^, and Di represents 

the first-derivative matrix with respect to the wall-normal direction (likewise for £>2 and 
£>3 below). The positive definite matrix Ad allows us to work in the £ 2 -norm. To do so, 
we apply a Cholesky decomposition to the weight matrix according to A4 = F'^ F and 
define ^ = £7 to arrive at 


£dV = I 'y^Ad'ydy= I F'^F-ydy = I = ||4||2 


( 2 . 21 ) 


where || • II 2 represents the £ 2 -norm and, accordingly, the eigenvalue problem (12.181) 
becomes 

- MF^) = FCF-\Fj). (2.22) 

Therefore, once the linear operator is redefined as Cl^ = FCF~^, we can conveniently 
use the £ 2 -norm and its associated inner product for all computations. The transient 
growth G, defined as the maximum energy growth over all possible initial conditions 
is given below in the £2 norm. 


G(t) = ||r ||2 = II 2 , (2.23) 

€0 || l ( 0)||2 |||( 0)||2 

where 7 ” is the linear evolution operator, i.e., the solution to equation (|2.16|) . 

The parameters that are to be investigated include the injection level G, the mobility 
parameter M, the charge diffusion coefficient Fe, the Taylor parameter T, the Reynolds 
number Re and the streamwise and spanwise wavenumbers a and j3. 


3. Numerical method and validation 

3.1. Numerical method 

To discretize the eigenvalue problem (I2.18L we use a spectral method based on colloca¬ 
tion points chosen as th e roots of Chebyshev polyn omials. The Matlab suite for partial 
differential equations by IWeideman fc Reddv ( 2000 f) is used for differentiation and integ¬ 
ration. 

To impose the boundary condition, we employ the boundary boarding technique ( Bovdl 


200 lh , in which selected rows of the linear matrices are replaced directly by the boundary 
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conditions. When solving the eigenvalue problem via the Matlab routine eig with the 
above boundary condition enforced, we find that the eigenvalues converge for a sufficient 
number N of collocation points (see figure IT^ and tabled in the validation section in the 
appendix!^ and approach the pure hydrodynamic results as electric effects become neg¬ 
ligible (see figure[T3]and table[6l). The corresponding eigenvectors, however, are incorrect 
since they do not satisfy the proper boundary conditions (not shown). To overcome this 
difficulty, we employ an iterative technique to obtain the eigenvector associated with a 
specified eigenvalue. In the generalized eigenvalue problem (12.161) , a desired eigenvalue uj 
(and its corresponding eigenvector) is targeted by applying the spectral transformation 


where 


S = (B - u;Ai)\Ai, 


—iA and S will be processed by an iterative routine (|Saadll201lll . 


(3.1) 


3.2. Validation 

The stability problem for EHD flow is exceedingly challenging from a numerical point of 
view, which warrants a careful and thorough validation step, before results about stability 
characteristics, modal and non-modal solutions and physical mechanisms are produced. 
To conserve the clarity of the paper structure, we postpone the validation steps in the 
appendix [A] 


4. Results of stability analysis 

4.1. EHD without cross-flow 

As mentioned earlier, the parameter T plays the main role of determining the flow in¬ 
stability. The critical Tc denotes the minimum value of T within the linear regime, above 
which infinitesimal disturbances can grow exponentially in time; will vary with the 
flow parameters. In the case of no cross-flow, the effects of Fe, T, M and C on the flow 
stability are investigated. As has already been assumed, the flow will be confined to the 
SCL (space-charged-limited) regime, implying a large value for C. 

We display the neutral stability curve in figure [3] for different Fe at C = 50, M = 100, 
T = 155, a = 2.5, and /? = 0. In the case without cross-flow, one does not need to 
distinguish between the x- and the z-axis, since neither is preferred by the base flow 
U = 0; thus, we simply set /3 = 0. As mentioned in the validation section, results for 
Fe = 10^ are very close to previous investigations. Even though the diffusion coefficient 
is small, it plays an important role in determining the critical T^, as shown in figure |3Da) 
and (b). For example, for Fe = 10^ the critical Tc decline s to 140. In fact, the valu e 
of Fe could fall within the range 10^ ^ 10^ for real liquids ( Perez fc Castellanoslll989l) . 
when Fe is nondimensionalized in the same way as presented here. Physically, the effect 
of diffusion will smooth out sharp gradients in the flow. Unlike the unidirectional electric 
field pointing in the wall-normal direction, the diffusion effect act equally in all directions. 
With charge diffusio n considered in the model, the discontinuous separatrix is blurred in 
the nonlinear phase ( Perez fc CastellanosI 1989ll . The physical mechanism of how charge 
diffusion influences the critical stability parameter Tc will be discussed by using an energy 
analysis (se e sectio nl5TJ) . In addition, the transient growth of disturbance energy has been 
discussed in AttenI ( I974I) using the quasi-stationary method; transient energy growth has 
been confirmed as a minor factor in this work. This is also confirmed in our computations, 
as presented in figure [3Kc): specifically, the figure shows that disturbance energy growth 
G reaches a value of about 3 at T = 155 for stable flows {Fe > 10^). 

The role of M in EHD is analogous to that of the Prandtl number in Rayleigh-Benard 
convection. In figure |4l(a), it is shown that the variation of M exerts no influence on the 
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Figure 3. Effect of Fe. The parameters are C = 50, M — 100, T = 155, a = 2.5, /? = 0 and 
N — 250. The direction of the arrow indicates increasing Fe. (a) Neutral stability curves for 
various Fe. (b) Tc as a function of Fe. (c) Transient energy growth versus time. 
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Figure 4. Effect of M. The parameters are C = 100, Fe = 10®, a = 2.57, /3 = 0 and N = 250 
without cross-flow, (a) The neutral stability curve, (b) Transient energy growth versus time for 
different M and T = 155. The direction of the arrow indicates increasing M. 


linear stability criterion, Tc = 159 .58 at C = 100, Fe = 10^, a = 2.57 and /3 = 0; the 
same finding has been reported in Atten fc Moreau ( 1972ll . For the transient dynamics, 
however, the same conclusion does not hold, as evidenced in figurelD^b). The plot describes 
a trend of increasing Gmax with smaller M. The slopes at the final time are slightly 
different for each M, indicating that the asymptotic growth rates differ slightly (while 
the linear stability criterion remains the same). 

Figure E] depicts the in fluence on C, which me asures the intensity of charge injection. 


Atten fc MorearJ ( 1972ll : Atten fc Lacroixl 1 1979ll reported a dependence of the critical 


value Tc on the parameter C. In figureEKa), we see that, in the SCL regime, increasing C 
will yield lower Tc. This result can be understood from a physical argument. Increasing 
the intensity of charge injection will lead to a higher concentration of charges between the 
electrodes. The linear instability mechanism, as discussed before, relies on the formation 
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Figure 5. Effect of C. The parameters are M = 100, Fe = 10"^, 13 = 0 and N = 250, without 
cross-flow. The direction of the arrow indicates increasing C. (a) Neutral stability curve for four 
different injection levels C within the SCL regime, (b) Transient energy growth for different C 
at r = 155, a = 2.57. 


of an electric torque due to convective motions. With higher charge concentration, the 
electric torque is stronger. Therefore, a lower voltage difference is required, which amounts 
to stating that a lower T will be sufficient to form an electric torque of comparable 
strength. But as we are in the SCL regime (with a value of C = 50 considered very 
large), a rise of C to 200 only yields a minor decrease in T^. In contrast, the transient 
dynamics of the perturbation energy G appears not to be influenced by a change in C 
for early times; for example, see the time interval t € [0, 3] in figure [SJd. 


4.2. EHD with cross-flow 

When cross-flow is considered, the property of the linearized system changes due to the 
presence of a base shear in the flow. Especially, this shear will render the linearized 
operator ‘more non-normal’. We first note that, in the modal stability analysis, Squire’s 
theorem still holds for EHD-Poiseuille flow, that is, a two-dimensional instability will be 
encount ered first. This can be easily verified by a perfect analogy with standard viscous 
theory ( Schmid fc Henningsonll200lll . Moreover, there are two sets of scales in the EHD 
problem with cross-flow. To study the influence of the cross-flow on the electric and the 
charge fields, the values of M and T are kept in the vicinity of the values in the previous 
section: the scale of the electric field will be considered primarily, whereas, when we 
examine effects of the electric field exerted on the canonical Poiseuille flow, we take 
the value of the free parameter Re around 5772, i.e., the linear stability criterion for 
pressure-driven flow; the latter choice introduces a scale based on the hydrodynamics. 
In both cases, we will enforce the relation Re = T/M'^, which results in the Reynolds 
number Re to be rather low in the former case (denoted as the low-i?e case) and relatively 
high in the latter case (referred to as the high-i?e case). 


4.2.1. EHD: low Re 

We have demonstrated that nonmodal effects in hydrostatic EHD are not significant. 
In the presence of cross-flow, given that the Reynolds number in this section is considered 
small, we expect the nonnormality to be rather moderate as well. For this reason, we will 
mainly focus on the modal stability characteristics for the low-Re case. 

In figure IHKa), it is observed that the symmetry of the hydrostatic EHD spectrum is 
now broken due to the presence of cross-flow. The most unstable perturbation travels at a 
positive phase speed Up = Wr/a = 2.25612.57 = 0.8778, induced by cross-flow convection 
(the centerline velocity of the cross-flow is 1, as we set Re = T/M"^ in equation (l2.9p L In 
figure EI d, we show the neutral stability curve for C = 50, M = 100, (3 = 0 and N = 250, 
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Figure 6. (a) Spectra for the case with and without cross-flow for C = 50, M — 100, Fe = 10®, 
T = 160, Re = T/M^ = 0.016 (for cross-flow), a = 2.57, /3 = 0 and N = 250. (b) Effect of Fe 
on the neutral stability curve. The parameters are C = 50, M = 100, /3 = 0 and N = 250, with 
cross-flow. The direction of the arrow indicates increasing Fe. 




Figure 7. Effect of M. (a) Neutral stability curve with cross-flow for C = 100, Fe = 10®, 
a = 2.57, 13 = 0 and N = 250. (b) Transient energy growth for different M and T = 145. The 
direction of the arrow indicates increasing M. 


which can be directly compared to the results in figure [3 Note that since Re = T/M^ is 
enforced, the Reynolds number Re is not identical for each point, but generally small. We 
see that, with cross-flow, the critical decreases compared to the no-cross-flow case; this 
indicates that the flow is more unstable in the presence of a low-Re cross-flow compared 
to the results in figure EKa) . To investigate the reason behind this destabilization, we 
again resort to an energy analysis in se ction 15.1.21 Previously^ an e nergy analysis for 
EHD with cross-flow has been studied in Castellanos fc Agrait ( 1992h . 

Even though varying M has no effect on the linear stability when C7 = 0, as has been 
discussed briefly in the previous section, in the presence of cross-flow, changing M does 
influence the linear stability. This is displayed in figure [T]) a), where we see that effects of 
M are only discernable when M is small. We will discuss this issue further in the energy 
analysis section 15.1.21 Considering non-normal linear stability, transient energy growth 
G is still small, even though slightly higher than in the no-cross-flow case. 


4.2.2. EHD: high Re 

In this section, we consider the flow governed by the inertial scale, i.e., in the high-i?e 
regime. To discuss the results more properly, the Reynolds number Re will be the free 
parameter, and the governing momentum equation is given by (12.101) . see section 1^751 The 
modal stability is examined in figure [8] In subfigure (a), changes in the spectrum due to 
the additional electric field are visible. It appears that the core modes, wall modes and 
center modes do not change appreciably, except that the growth rate of the most unstable 
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Re 


Figure 8. (a) Comparison of two spectra for T = 10 ® and T = 200 at 
C — 100,-Re = 5600, Fe = 10®, a = 1,/1 = 0. (b) Effect of T on the neutral stability curve 
with cross-flow at C = 100, Fe = 10®, M — /3 = 0. The direction of the ar¬ 

row indicates increasing T. (c) Effect of Fe on the neutral stability curve with cross-flow at 
C = 100, T = 100, M = ^JT /M^, /3 = 0. The direction of the arrow indicates increasing Fe. 


mode increases. In subfigure (b), we plot the neutral stability curve for varying T. The 
pure hydrodynamic linear stability limit Re = 5772.2 is recovered by considering a minute 
value for T such as T = 10“® (we could have taken T = 0, but to be compatible with 
equation (12.81) and the discussion based on that equation in other literature, we assign 
to T a negligibly small value). With increasing T, the system becomes more unstable, as 
the critical linear stability criterion becomes smaller. The reason for this is obviously due 
to the effect of the electric field transferring energy into the velocity fluctuations, while 
at the same time modifying the canonical channel flow; see table [3] in the energy analysis 
section IS. 1.31 at C = 100, .Fe = 10®, .Re = 5500, a = 1 and /3 = 0. 

We also investigate the effect of charge diffusion Fe on the flow stability. The results 
concerning the neutral stability curve are shown in figure[8l(c) at C = 100, T = 100, M = 
13 = 0. For small charge diffusion (large Fe), the critical Reynolds number Re is 
only slightly affected by changes in Fe. Only when Fe = 10® does the critical Reynolds 
number Re drop noticeably, though the destabilization effect is still small. It thus can 
be concluded that charge diffusion has only a small influence on the dynamics of EHD 
cross-flow at high Re. This is due to the inertial scale we are considering. As we have 
seen in the hydrostatic EHD flow, the effect of charge diffusion is significant, considering 
the electric scale, i.e., at relatively small (or zero) Reynolds numbers. 

It is well known that in high-.Re Poiseuille flow the two-dimensional Orr mechanism is 
not the principal mechanism for perturbation energy growth over a finite time horizon. 
The flow is expected to become turbulent within a short time interval, even though the 
asymptotic growth rate of the linear system is negative. The non-normal nature of the 
linearized Navier-Stokes operator for channel flow — in physical terms, due to the base 
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Figure 9. Effect of T. The direction of the arrow indicates increasing T. (a) Transient energy 
growth for C = 100, Re = 5200, Fe = 10®, a — 0, S = 2. (b) Contours of transient growth G in 
the a-/l-plane at C = 100, Re = 5000, Fe = 10®, T = 100. (c) Transient growth G as a function 
of /3 at G = 100, Fe = 10®, a = 0. 


flow modulation by spanwise vorticity tilting into the wall-normal direction — snggests 
that transient disturbance growth during the early phase should be considered primarily. 

In figure |9l[a) , we present the transient growth G for different T. Mainly, the effect of 
increasing T is to enhance transient growth. The optimal initial condition which achieves 
maximum transient growth is shown in figurelH^b). The optimal wavenumbers for the pure 
hydrodynamic case, independent of Re, are found to be a = 0 a nd 0 = 2.05, suggesting 
streamwise-independent vortices as the most amplified structures (jSchmid fc Henningson 
l200lh . For high-i?e EHD with cross-flow, the maximum transient growth is still found 
to favor streaks (a = 0), but with a different optimal spanwise wavenumber of /3 = 2.36 
at C = 100, i?e = 5000, Fe = 10®, T = 100 (see figure [^Kb)). Interestingly, for a different 
value of T, i.e., a different amount of potential drop across the electrodes, the optimal 
wavenumber would be different. For example, at T = 50 the optimal /3 = 2.18, while 
at F = 10 the optimal (3 = 2.08, as shown in figure [91(c) . It seems that, for smaller 
T, approaching the regime of pure hydrodynamics, the optimal P is converging towards 
P = 2.05. The independence of maximum transient growth on Re for various P still holds 
in the limit of high-i?e EHD flow, as indicated by the dashed lines connecting the peaks 
of the two curves for Re = 5000 and Re = 3000 in figure (HKc). In the nonlinear regime of 
E HD Poiseuille flow, the infl uence of the electric field on the streaks has been reported 
in Soldati fc Baneriee ( 1998ll . These authors reported the spanwise spacing of the low- 


speed streaks are about 105±15 in wall units, which is different from the average spanwise 
spacin g of the streaks in Poiseuille flow, that is 100 in wall units, see iButler fc Farrell 
1 19931) for example. Thus, to some extent, our results that the spanwise spacing of the 
streaks changes in the linear EHD cross-flow agree qualitatively with these findings. 
These authors also found that the cross-flow is weakened by the electric field. This does 
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Figure 10. Effect of T. (a) The optimal initial condition for t) as a function of T and (3 at 
C = 100, Fe — 10®, Re = 5000, a = 0. (b) The optimal initial conditions of velocity and 
potential at C = 100, Re = 5000, Fe = 10®, T = 100, a = 0, = 2.36. (c) Velocity vectors of 

the optimal initial condition in the cross-stream plane; same parameters as in (b). (d) Contour 
of the optimal response for (/?; same parameters as in (b). 


not stand in conflict with the current results, because enhanced transient growth due 
to the electric field, as found here, only indicates that, in the linear phase, transition 
to turbulence is more rapid when compared to canonical channel flow; no conclusions 
can be drawn for the flow behavior in the nonlinear regime. To more fully understand 
how the electric field influences streaks and streamwise vortices in the nonlinear phase 
of transition, a more comprehensive st udy of the role played by EHD in the formation 
and dynamics of a self-sustaining cycle ( Jimenez fc Pinellilll999l l is called for. 

To further investigate the effect of T, we plot in figure [TO] the optimal initial conditions 
which achieve Gmax in a given finite time for parameters C = 100, Re = 5000, Fe = 10®, 
a = 0 and different values of T with its corresponding optimal /3. In subfigure (a), 
the optimal initial conditions for v for various T are presented. The symmetry of the 
optimal V with respect to the flow centerline y = 0, when the flow is close to the pure- 
hydrodynamics limit, is broken due to the action of the electric field in the wall-normal 
direction as T increases. Since the electrode with higher potential is at ?/ = 1 in our case, 
the optimal initial conditions for w and v are tilted towards y = 1 (see also subfigure (b), 
which additionally shows the optimal initial condition for ip). In subplot (c), the formation 
of streamwise vortices in the y-z-plane is shown; their centers are shifted upwards by the 
electric field. In subfigure (d), the optimal response of ip^ taking the form of waves in the 
spanwise direction, is displayed. Recalling that the nonmodal transient growth is due to 
base-flow modulations arising from the tilting of spanwise into wall-normal vorticity, we 
can state that the variation of the optimal spanwise wavenumber for different T is the 
direct result of the three-dimensional nature of the non-normal linearized operator under 
the influence of a constant electric field pointing in the wall-normal coordinate direction. 
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5. Results of energy analysis 

5.1. Asymptotic energy analysis 

The dynamics of the disturbance energy (of the velocity fluctuations) in the limit of 
an infinite time horizon is examined in this section. The governing equation for the 
energy evolution is obtained by multiplying the linearized equation (I2.8b l by the complex 
conjugate velocity vj, i.e., 

dtp ^ 
dxjXj dxi'’ 

(5.1) 

taking the complex conjugate of the obtained equation, and averaging the two equations, 
which leaves us with 


t I t 


dxi 


trr 


- 


dp 

dvi 




d^Vj 


T 


dxj 
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2 dxiXi 


(5.2) 


where £ = v\vil‘l is the perturbation energy density of the hydrodynamic part in the 
spectral space. The terms in the square brackets are the transport terms which, in case 
of periodic as well as no-slip and no-penetration boundary conditions , exert no influence 
on the energy balance. Therefore, after integrating the above equation over the control 
volume O, we obtain 


dt 


dV = - 


R t ^ t .w f dv. 
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(5.3) 


VE3 


Since the boundary conditions are periodic in the wall-parallel coordinate directions, it is 
legitimate to consider the control “volume” fl only in the j/-direction, that is, fl = [—1,1] 
and dV = dy. The first term on the right-hand side of equation (15.31) represents energy 
production from the mean shear (Pr), which is zero in the hydrostatic case; the second 
term describes viscous dissipation (VD); the third to hfth terms are the energy transfer 
terms between the velocity fluctuation field and the electric field (VEl, VE2, VE3, respect- 
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ively). The Einstein summation convention does not apply for the subscripts of VElij, 

, , r ■ I “-r dip''^ dv d(p dv^ 

example; the term VEI 21 , tor instance, represents ' 


dcj) 


(9.x dx dx (9.x 


for 
dV. 


As has been discussed and verified for polymeric flows in Zhang et al. ( 2013ll . the time 
variation of the normalized perturbation energy density should be equal to the twice the 
asymptotic growth rate of linear disturbances, i.e.. 


Re = 


5 - 


- 2 (^ 2 , 


(5.4) 


£dV 


where denotes the growth rate of the least stable mode. We will validate this relation 
in the following sections and use it as an a posteriori check for our results. 


5.1.1. EHD without cross-flow 

We apply the energy analysis for hydrostatic EHD flow with different values of charge 
diffusion coefficients Fe to probe how the electric field interacts with the velocity fluctu¬ 
ations. Quantitative results are listed in table [U with the notation VD = VDn -t VD12 -I- 
VD21 + VD22 and, likewise, VE = VEI21 + VEI22 + VE2 +VE3. No spanwise dependence as 
/3 = 0. Immediately, one can make several direct observations. First, viscous dissipation is 
always negative for the hydrodynamics. Second, since the EHD flow is hydrostatic, there 
is no production from the mean shear, Pr = 0. The only terms that can lead to growths in 
the hydrodynamic disturbance energy density £ are linked to the energy transfer from the 
electric field, VE. The most efficient mechanism seems to be related to the term VEI21, 
which represents the interaction between the streamwise perturbed electric field and the 
wall-normal velocity shear under the constant effect of the wall-normal base electric held. 
The term VE3 is even negative, indicating that the electric held can absorb energy from 
the perturbed hydrodynamic held by an out-of-phase conhguration between ip and v (in 
the energy-budget equation for the perturbed electric held, one would hnd the exact 
same term with opposite sign). Regarding the effect of charge diffusion, with increasing 
1 /Fe (increasing charge diffusion) from the right column to the left in the table, the total 
energy transfer VE diminishes, but, at the same time, the hydrodynamic diffusion is also 
dissipating less energy into heat. Furthermore, even though VE (and thus VEI22, VE2 and 
VE3) decreases with rising charge diffusion, the primary mechanism of energy transfer 
VEI21 transfers more energy from the electric held to the hydrodynamic huctuations, 
together with a less dissipation, leading to an unstable how for the chosen parameters. 
Therefore, it seems that the effect of charge diffusion is to catalytically enhance the 
efhciency of the most productive energy transfer mechanism between the perturbed elec¬ 
tric held and the hydrodynamics, expressed by the term VEI21, and result in a lower 
dissipation. As a consequence, increasing charge diffusion leads to a more unstable how. 

It is instructive to assess the effect of M on the linear stability (see section HIT]) with the 
help of the energy-budget equation (15.31) . For a vanishing time derivative of the energy 
density, the factor on the right-hand side can be eliminated for the hydrostatic case 
U = 0. In the case of cross-how, however, with the same reasoning M will have an 
inhuence on the linear stability criterion. 


5.1.2. EHD with low-Re cross-flow 

Table [2] shows the results for C = 50, M = 100, T = 160, Re = T/M"^ = 0.016, 
a = 2.57, 13 = 0 with cross-how for different Fe. Compared to the case without cross-how, 
the results are quite similar. However, it is interesting to note that the huctuation energy 
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terms 

Fe = 10^ 

Fe = 10® 

Fe = 10® 

Fe = 10 '^ 

VDii 

-254.9836 

-255.3751 

-255.4628 

-255.4781 

VD 12 

-481.9379 

-486.1119 

-487.0231 

-487.1773 

VD 21 

-570.6289 

-570.2374 

-570.1497 

-570.1344 

VD 22 

-254.9836 

-255.3751 

-255.4628 

-255.4781 

VEI 21 

1154.1206 

1152.8793 

1152.6403 

1152.5936 

VEI 22 

428.695 

431.8134 

432.4984 

432.6197 

VE2 

30.0141 

30.7628 

30.9299 

30.9614 

VE3 

-50.2196 

-48.3506 

-47.9797 

-47.9186 

VD 

-1562.5341 

-1567.0994 

-1568.0984 

-1568.2679 

VE 

1562.6102 

1567.1049 

1568.089 

1568.2561 

Pr 

0 

0 

0 

0 

Re 

0.0761 

0.0054 

-0.0094 

-0.0119 


0.0761 

0.0054 

-0.0094 

-0.0119 


Table 1. Energy budget for modal instability of hydrostatic EHD flow for different values of the 
charge diffusion. The results are normalized as in equation dEl. The parameters are C = 50, 
M = 100, T = 160, a = 2.57, /3 = 0 for hydrostatic flow. 


terms 

Fe = lO"* 

Fe = 10® 

Fe = 10® 

Fe = 10'^ 

VDii 

-254.1325 

-254.1337 

-254.1318 

-254.1309 

VD 12 

-475.0779 

-476.292 

-476.5068 

-476.5348 

VD 21 

-571.48 

-571.4788 

-571.4807 

-571.4816 

VD 22 

-254.1325 

-254.1337 

-254.1318 

-254.1309 

VEI 21 

1158.7772 

1157.62 

1157.3579 

1157.3025 

VEI 22 

422.786 

423.6012 

423.7687 

423.7989 

VE2 

26.3245 

26.363 

26.3747 

26.3794 

VE3 

-52.8799 

-51.3981 

-51.1094 

-51.0629 

VD 

-1554.8229 

-1556.0382 

-1556.2511 

-1556.2782 

VE 

1555.0078 

1556.1861 

1556.3919 

1556.4179 

Pr 

-0.0079038 

-0.0080187 

-0.0080507 

-0.0080566 

Re 

0.1769 

0.1399 

0.1327 

0.1316 

2u]i 

0.1769 

0.1399 

0.1327 

0.1316 


Table 2. Energy budget for modal instability of EHD flow with low-7?e cross-flow for different 
values of the charge diffusion Fe. The results are normalized as in equation (15.411 . The parameters 
are C = 50, M = 100, T = 160, Re — T/M^ = 0.016, a = 2.57, /? = 0 with cross-flow. 


production from the mean shear Pr, even though rather small, is negative, indicating that 
the perturbed flow field transfers energy to the base flow. Recalling the results in figure 0 
of section 14.2.11 a change of M does not have a strong effect on the rate of change of the 
disturbance energy density E since Pr is very small. 

In the case of EHD flow with a weak cross-flow {Re = 0.016), the main mechanism 
for transferring energy into the hydrodynamic subsystem is still based on the potential 
difference across the two electrodes — the same as for the no cross-flow case. This can 
be confirmed by inspecting tabled VEI 21 is the dominant energy transfer term. 

5.1.3. EFID with high-Re cross-flow 

The energy analysis for the EHD Poiseille flow at C = 100, Fe = 10^, Re = 5500, 
a = 1 and /3 = 0 is summarized in table[3l Note that the production Pr is dimishing with 
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terms (xlO 

T = 10"® 

T = 10 

T = 100 

T = 200 

VDii 

-2.643 

-2.643 

-2.6429 

-2.6429 

VD 12 

-138.22 

-138.24 

-138.42 

-138.61 

VD 21 

-0.99338 

-0.99339 

-0.99343 

-0.99347 

VD 22 

-2.643 

-2.643 

-2.6429 

-2.6429 

VEI 21 

3.1782 

0.31785 

3.1814 

6.3689 

VEI 22 

7.3458-IQ-^ 

0.73467 

7.355 

14.729 

VE2 

7.0501 -lO-^ 

0.070507 

0.70555 

1.4122 

VE3 

-1.2043-lO"^® 

-0.12044 

-1.2055 

-2.4133 

VD 

-144.5 

-144.52 

-144.7 

-144.89 

VE 

10.025 -lO"^® 

1.0026 

10.036 

20.096 

Pr 

132.35 

132.06 

129.41 

126.44 

R. 

-12.153 

-11.463 

-5.2526 

1.6542 


-12.153 

-11.463 

-5.2526 

1.6542 


Table 3. Energy budget for modal instability of EHD flow with a high-_Re cross-flow for dif¬ 
ferent values of the stability parameter T. The results are normalized as in equation (15.411 . The 
parameters are C = 100, Fe = 10®, Re = 5500, a = 1 and (3 = 0 with cross-flow. 


increasing T. On the other hand, VE increases with larger values of T, compensating 
and exceeding the decrease of Pr at higher T. This is consistent with the results in 
figure[5] higher values of T yield a more unstable flow. However, the principal mechanism 
underlying the flow instability is still linked to the production Pr. The electric field only 
assumes a secondary role in destabilizing the flow, at least for the parameters considered 
in this case. Unlike the hydrostatic case where VEI 21 is responsible for the dominant 
energy transfer, in the presence of cross-flow VEI 22 becomes the most efficient agent 
transferring fluctuation energy £ between the electric field and the perturbed velocity 
field. 


5.2. Transient energy analysis 

To investigate the cause for the increase of nonmodal growth with T, recall figure [H 
we formulate and perform an energy analysis for the initial-value problem of equa¬ 
tion 112.161). We consider the energy density evolution over a finite time horizon fol¬ 
lowing Butler &: Farrelll ( 1992ll 


1 

M 


f d£ 1 p d f \ 

~W\J-Jo Jo 2 ) 


dV 


1 1“^ ' dU 1 dui dui 

MJ~iJo Jo V dy Redxjdxj 


Pr VD 


T d(j) dip dv T d^(j) dp 
Re dy dxj dxj Re dy'^ ^ dy 

VEl VE2 


(5.5) 


T d^4) ' 


VE3 


dz dx dy 


where |H| = 2ab, a = 2 -k ja and b = 21:jf. In the above equation, we label, as before, the 
first term on the right-hand side as Pr (production from the mean shear), the second term 
as VD (viscous dissipation), the third to fifth terms collectively as VE (energy density 
transfer between the perturbed velocity held and the perturbed electric held) and the 
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Figure 11. Energy analysis over a finite time horizon at C = 100, Re — 5000, Fe= 10®, 
a = 0. (a) T = 10“®, /3 = 2.05; (b) T = 100, /3 = 2.36. 


sum of all five terms as Total. In this temporal evolution problem, the initial condition 
is the optimal one, following the procedure in section 14.2.21 

Results of our energy analysis are presented in hgure 1111 In subfigure (a) and its 
inset, the pure hydrodynamic result is shown, where the production Pr counteracts the 
viscous dissipation VD. In subplot (b) for EHD cross-flow, we observe that the term VE is 
insignificant, even though at T = 100; this is in contrast to both the linear modal stability 
criterion (see figure [S]) and the overall nonmodal transient growth (see figures IH] and ITUl) 
where its effect is not negligible. Furthermore, production Pr increases by a factor of 2 ~ 3 
compared to the pure hydrodynamic flow. These results seem to indicate that, concerning 
the nonmodal analysis, the effect of the additional electric held on the canonical channel 
flow is incidental, i.e., the perturbation velocity energy is only indirectly inhuenced by 
the electric held; in fact, the electric held does not induce a substantial energy transfer 
directly into velocity huctuations at all and its effect is to enhance the lift-up mechanism, 
therefore the production. 

Examining more closely the inset in hgure ITlT bl. we see that the term VE surpasses Pr 
only in the very beginning of the time horizon. This is due to the high-i?e regime we are 
investigating. As discussed earlier. Re = T /represents the ratio of the momentum 
relaxation time to the charge relaxation time L*'^/{K*A(j)Q). With the maximum 

Total energy achieved at tmax ~ 144 in hgure fTTI bl. we can estimate the time horizon in 
the inset by observing that 144/5000 = 0.029, a value close to the time scale depicted in 
the inset of (b). Moreover, the minimal energy growth due to the electric force validates 
our previous observation that the purely EHD-induced non-normality is rather small (see 
appendix |B] for a direct proof of this statement via an input-output formulation). In the 
case of other complex hows at high Reynolds number s , a similar con clusion can be draw, 
for instance, in viscoelastic hows ( Zhang et al. 2013 : Brandt 2014l l. polymer stretching 
cannot induce disturbance growth when the huid inertia is prevalent. 


6. Discussion and conclusions 

In this article, we have presented a comprehensive linear stability analysis of charge- 
injection-induced electrohydrodynamic hows between two plate electrodes, covering the 
hydrostatic as well as the cross-how case, employing modal as well as nonmodal tools. We 
intend to examine whether the linear framework is sufficient for describing the transition 
to turbulence of EHD how in the early phase of perturbation evolution. It is hoped that 
the results presented above and summarized below would help to understand better the 
EHD how instability and its transition mechanism and shed light on its how physics as 
well as how control design. 
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6.1. EHD without cross-flow 

In the hydrostatic case, the often-omitted charge diffusion is taken into account and found 
to have a non-negligible effect, particularly on the critical linear stability parameter 
with SCL injection — a finding running contrary to a common assumption in previous 
studies. In those studies, a linear stability analysis predicts a critical value of Ri 161 in 
the strong injection limit. This result is reproduced in our computations for a negligible 
value of IfFe, but even for a moderate amount of charge diffusion the flow quickly 
becomes more unstable. Hence, we suggest that charge diffusion be accounted for in 
linear stability analyses and numerical simulations whenever the real physics indicates 
charge diffusion that can not be neglected, as it improves not only the model of the flow 
physics but also the robustness of the numerics as we ll. In fact, the common use of total 
variation diminishing (TVD) schemes (lHartenlll983^ in direct numerical simulations of 
EHD flow, which introduces artificial numerical diffusion, seems unnecessary when true 
physical charge diffusion could be included. 

The longstanding discrepancy of the critical stability parameter between the ex¬ 
perimental and theoretical value, however, could not be resolved by our analysis: even 
th ough Tr, in the SCL limit drops to 140 at Fe = 10^ (a physical value according 


to 


Perez fc Castellanos! ( 1989lH . a substantial gap remains to the experimentally meas¬ 
ured parameter of Tc « 110. Motivated by the researches in subcritical channel flow, 
we examine other mechanisms for early transition to turbulence, specifically, transient 
growth due to the non-normality of the linearized EHD operator. Nonmodal stability 
theory has been successfully applied to the variety of wall-bounded shear flows in an 
attempt to explaining aspects of the transition process. In the case of hydrostatic EHD 
flow, our calculations seem to indicate that transient energy growth, as defined in equa¬ 
tion (j2.19|) . is not significant, reaching gains of ^ 10 at most: the flow instability is rather 
dictated by the asymptotic growth rate of the least stable mode. These results seem to 
indicate that below the critical Tc the significant energy growth observed in the real EHD 
flow is not of a linear nature, otherwise the linear framework would succeed to detect 
it. It might be hypothesized that the major energy growth mechanism in subcritical hy¬ 
drostatic EHD follows a nonlinear route; nevertheless, it is only after performing a full 
nonlinear simulation of subcritical EHD flow that can one conclude whether its energy 
growth mechanism is truly nonlinear or not. Besides, these results also seem to shed some 
light on the flow control of hydrostatic EHD. It is now well-established that in canon¬ 
ical channel flow, where the perturbation energy is found to grow linearly in the early 
phase, a linear flow control stra tegy is sufficient to abate the perturbation development 
(lKimll2nn3l: rK im fc Bewlevll2007h . Due to the limited early perturbation energy growth in 


hydrostatic EHD, we thus suggest that different flow control methods be examined and 
applied in addressing the flow control of such flo w. There migh t exist another possibility 
for the limited transient growth. As discussed bv I Atteii (1974), the correct prediction of 
the linear stability criterion might require a closer comparison between the experimental 
conditions and the mathematical model. In this light, one may suggest a re-examination 
of the charge creation and transport processes, as the current charge creation model does 
not seem to accommodate any efficient energy transfer from the electric to the flow field, 
during the linear phase. 


6.2. FED with cross-flow 

The flow instability and the transition to turbulence in canonical or complex channel (for 
example, EHD, MHD or polymeric) flows are currently not well understood. The study 
on the complex channel flow, serving as a supplement to the investigation of the canonical 
flows, focuses on the flow modification under the influence of external fields, for example. 























24 


M. Zhang, F. Martinelli, J. Wu, P.J. Schmid and M. Quadrio 


electric field, magnetic field or polymer stress field. The study of such flows will not 
only improve our understanding of these particular flow configurations, but also, more 
importantly, help us to better understand, during the linear, transition and turbulent 
phases, the dynamics of the important flow structures, for instance, the streak formation 
and attenuation, by probing the interaction between the fluids (or flow structures) and 
the external fields. For example, in the point of view of flow control, the research on 
polymer turbulence drag reduction reveals the mechanism how t he auto-generation cycles 
of turbulence is modified in the presence of polymer molecules ( Dubief et al.ll2003l . This 
has led to an even broader picture of the dynamics of turbulence. Similarly, in the case 
of EHD, our goal is to understand how the flow changes in response to the electric effects 
and provide a physical interpretation. Below we present the results of EHD cross-flow. 
We differentiated low-i?e and high-i?e cases. For low-iie flow, the effects of M and Fe are 
similar to those of the hydrostatic flow, with the linear stability criterion being smaller 
for low-i?e cross-flow when compared to hydrostatic flow. 

The high-i?e case is more interesting. In both modal and nonmodal stability ana¬ 
lyses, the canonical channel flow becomes more unstable, once an electric field is applied 
between the two electrodes. From an input-output and an energy analysis we found, 
however, that the energy growth directly related to the electric field is not significant 
and that the effect of the electric field on the flow instability is indirect. In general, in 
high-i?e channel flow, the maximum transient growth is achieved by vortices aligning 
along the streamwise coordinate direction and generating streamwise streaks via an ef¬ 
ficient energy growth mechanism known as lift-up. These optimal streamwise vortices 
are symmetric with respect to the channel centerline for standard Poiseille flow. In con¬ 
trast to other complex flows, in EHD flows the electric field, which always points in the 
wall-normal direction, actively participates in the formation of the streamwise rolls by ac¬ 
celerating the downward-moving fluid (note that, in our setting, the injector is at ?/ = 1). 
Consequently, this yields stronger transient growth via the lift-up mechanism, when com¬ 
pared to the common channel flow. In other words, the electric field provides wall -normal 
mome ntum. As has been discussed in iLandrihl I I980l l and recently reviewed by [Brandt 
( 2 OI 4 I I. the presence of wall-normal momentum will cause any three-dimensional, asymp¬ 
totically stable or unstable shear flow to exhibit energy growth during a transient phase. 
In the present study, the role of the electric field is to provide the shear flow with such 
a source of wall-normal momentum and to strengthen the lift-up mechanism for EHD 
flow with high-i?e cross-flow. Besides, we also find that the optimal wavenumbers for 
maximum transient growth increase under a stronger electric effect. Since the electric 
field will help to establish streamwise vortices, it may constitute a good actuator for 
drag reduction techniqu es, using the two-dimensio n al rolls together with a flo w control 
strategy as described in Schoppa fc Hussain ( I998ll ; Soldati fc Baneri^ ( 1998 ). 
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Appendix A. Code validation 

We first perform a resolution check to examine the convergence of the results. The 
parameters in this case are C = 50, Fe = 2000, Re = 6000, T = 100, M = ^JT/Re = 
0.129, a = 1 and /? = 0. The eigenspectra for four different grid resolutions N are shown 
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Figure 12. Code validation, (a) Resolution check for EHD flow with cross-flow at C = 50, 
Fe = 2000, Re = 6000, T = 100, M = ^T/Re = 0.1291, a = 1 and f = 0. (b) Eigenvector 
component v for the most unstable mode (normalized to have the same peak value for the two 
codes). The parameters are C = 50, M = 10“^^, Fe = 2000, Re = 6000, T = ■ Re, a = 1, 

P = 0 and N = 250 for the EHD code, and Re = 6000, ot = 1, f = Q and N — 250 for the 
hydrodynamic stability code. 


N most unstable mode 


150 0.260023637950300 -f 0.000652797815269i 

200 0.260023637089882 -f 0.000652796819289i 

250 0.260023637069851 -f 0.000652796791624i 

280 0.260023637052960 -f 0.000652796810920i 

Table 4. Code validation. Resolution check for the most unstable eigenvalue of EHD flow 
with cross-flow; with same parameters as in figure fT^ 


in figure [T^ a). The most unstable mode in these cases are listed in table ID Satisfactory 
convergence, with increasing N, is observed. 

Secondly, the EHD eigenvector, from using is examined against a verified, pure 

hydrodynamic stability code employing the same spectral collo c ation method and solving 
the Orr-Sommerfeld-Squire system, see lSchmid fc Henningson ( 2001 1. as shown in figure 
[T2](b). The parameters for the EHD code are C = 50, M = 10“^^, Fe = 2000, Re = 6000, 
T = • Re, a = I, (3 = 0 and N = 250 . The parameters for the hydrodynamic 

stability code are Re = 6000, a = 1, f = 0 and N = 250. We see that the iteratively 
solved EHD eigenvector is the same as the pure hydrodynamic one, which is solved 
directly by the Matlab routine eig. For the computation of the transient amplification 
G in e quation (12.231). it is legitima te to include only the first several, most unstable 
modes (|Schmid fc Henningsonll2001 1. i.e., eigenmodes corresponding to eigenvalues with 
imaginary part smaller than a certain wf are discarded, see table [5] for a validation of 
this approach. The reason for a minor increase of G, as more modes are included, is due 
to the newly incorporated eigenvectors, not because of an insufficiently refined grid. 

With the eigenvalue problem reliably solved as shown above, we present validation for 
the specific flows considered here. In the case of hydrostatic flow, our results for Fe = 10^, 
approximating the case of zero charge diffusion, = 160.67 and a = 2.57 at C = 50 
(see figure Ha) in section 14.11) . are very close to the linear stability criterion reported 
in Atten fc Moreau ( 1972h . Tc = 160.75 and a = 2.569 in the case of C —>■ oo, where 


a coupled flow and electric system with neglected charge diffusion has been considered. 
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q G ujt G 


-3.5 3.312404 -0.5 1.138543e + 04 

-5.5 3.434333 -3.5 1.174459e + 04 

-7.5 3.434351 -10.5 1.175565e + 04 

Table 5. G versus different cut-off growth rates ojf. The first two columns represent hydrostatic 
EHD flow at C = 50, Fe = 10®, M = 100, T = 155, a = 2.5, /3 = 0, = 250. The last two 
columns represent EHD flow with cross-flow at C = 100, Fe = 10®, Re = 5000, T = 100, a = 0, 
P = 2.36, N = 250. Refer to the text for the definition of cuf. 


EHD code hydrodynamic stability code 


^max,r 

0.259815871017297 

0.259815871062631 

^max,i 

0.000323088678313 

0.000323088655527 

^max 

18.86745 

18.87514 

Gmax 

38.92401 

38.93307 


Table 6. Code validation. The parameters are the same as in figure [TS] 


This additionally implies that a value of C higher than 50 can well approximate the 
space-charge-limit. 

In the presence of cross-flow, since there exist no quantitative results for eigenvalues 
and eigenvectors of the EHD problem in the literature, we partially verify our results 
by examining the pure hydrodynamic limit of the EHD-linearized problem, i.e., with 
electric effects being very small. This comparison is made with the stability code. The 
parameters are identical to the ones chosen above for the comparison of the eigenvectors. 
The Poiseuille base flow \s U = 1 — in both codes. It is obvious that, with these 
selected parameters, the governing equation (12.13b ') for Lp is void of the coupling with 
V, since in equation (I2.131 i,i is negligible. Therefore, the hydrodynamics equations 
for V and rj in (12.13b ') and (12.13b ') must reproduce the results of the stability code. 
This match is shown in figure [121 In subfigure (a), the spectra of two codes are seen to 
collapse, even in the intersection region of the three eigenbranche s, which is known to be 
sensi tive due to the high non-normality of the linearized system ( Schmid fc HenningsonI 
200lh . Additionally, the blue eigenmodes in (a) not matched by the red hydrodynamic 


modes are the supplementary eigenvalues linked to the presence of an electric field. The 
most unstable eigenvalue is shown in table [51 In subfigure (b), transient growth using an 
eigenvector expansion with n = 71 eigenmodes is shown. A quantitative comparison of 
the maximum transient growth G^ax ^^nd its corresponding time t^ax is presented in in 
tableini Agreement up to the fourth digit is achieved. The computations of tmax and G^ax 
involve n = 71 eigenfunctions, each one solved with the iterative method. Even though 
each individual mode may be prone to small inaccuracies, figure I13l bi illustrates that 
transient growth (a multi-modal phenomenon) can be reliably and robustly computed 
using the eigenvector expansion outlined above. 


Appendix B. Input-output formulation 

An input-output formulation can reveal additional information on prevalent instabil¬ 
ity mechanisms by considering different types of forcings (input) and responses (out- 
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Figure 13. Comparison between EHD-Poiseuille flow and canonical Poiseuille flow. The para¬ 
meters are C = 50, M = 10“^^, Fe = 2000, Re = 6000, T = - Re, a = 1, = 0 and N = 250 

for the EHD code, and Re — 6000, a = 1, /3 = 0 and N — 250 for the purely hydrodynamic 
stability code, (a) The eigenvalue spectrum, (b) Transient amplification of initial energy. 


input Bi„i{v,r],ip} Bin 2 {v,y} Bi„i{ip} 


583.4087 587.7147 552.0443 

11765.40 11350.38 428.8926 

Table 7. Result from an input-output analysis at C = 100, Fe = 10®, Re = 5000, T = 100, 

a = 0 and /3 = 2.36. 


put) ( .Tovanovic & BamiebI 2005 1. To demonstrate that the transient growth due to per¬ 
turbative Lp is small, we compare the full responses to perturbations consisting of (i) all 
variables v,rj and (/?, (ii) both v and rj, and (hi) only p. We thus dehne for these three 
cases different input filters J5, where Bmi = Isnxsn for the first case, while for the 
second and third cases we have 



(Bl) 


The output filter Cout is Isnxsn for all three cases: we examine the flow response in 
all velocities and the electric field. Consequently, the energy weight matrices should be 
redefined with Adout = CoutMCj^t and Mi„ = After applying a Cholesky 

decomposition to these energy weight matrices, we obtain Font and Fin for a formulation 
based on the L 2 -iioriR- Finally, the maximum transient growth G over a finite time 
interval is given by 


G{t) = max 


||'7o'ut(0l|So 


"I'D ll7»n(0)llK 

FoutTF-^ 


= max 


Iir7« 
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= max- 
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= max- 
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\Fin'^ii 


\Fin'^ii 


= WFoutTF, 


1||2 
in \\^ 
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tc. 


BinF^^ II2. 
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We report the transient growth results for the above three cases in table [7] at C = 100, 
Fe = 10®, Re = 5000, T = 100, a = 0 and /3 = 2.36. We observe that perturbations 
solely in ip (case (iii)) exhibit transient growth two orders smaller than in the other 
two cases. For cases (i) and (ii) the transient growth characteristics are nearly identical 
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which suggests that the nonnormality of the linear operator is mainly related to the 
hydrodynamics. 
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